clear;
clc;
close all;

%% 参数设置
Fs = 96e6;              % 采样率96MHz
f_signal = 400e3;      % 信号频率400KHz
t = 0:1/Fs:1e-3;       % 生成1ms时长信号（96000点）
s = sin(2*pi*f_signal*t);  % 原始正弦波
x = awgn(s, 40);        % 添加SNR=30dB的高斯噪声
D = 24;         % 降采样因子
N = 5;          % 5阶CIC
input_bit = 4;  % 输入4bit无符号数
output_bit = 4+ceil(N*log2(D));

%% SDADC滤波器处理
num = [1, -2.994219815030523, 2.994219815030523, -1];
den = [1, -2.33063271405379, 1.866676648819709, -0.508718943549610];
y = filter(num, den, x); % 滤波得到SDM输出

%% 频谱分析
Y = fft(y);
P = abs(Y)/length(Y);
f = Fs*(0:length(Y)-1)/length(Y);
% 绘制单边频谱（0 ~ Fs/2）
N_half = floor(length(Y)/2);         % 取前一半频谱
figure;
plot(f(1:N_half), 20*log10(P(1:N_half)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('SDADC输出频谱（单边）');
grid on;
xlim([0, Fs/2]);

%% 量化到4bit无符号数
% 将y归一化到[0, 15]范围
y_normalized = (y - min(y)) / (max(y) - min(y)); % 归一化到0~1
y_quantized = round(y_normalized * 15);       % 映射到0~15（4bit无符号）
% y_quantized = uint8(round(y_normalized * 15));   % 4bit量化
% y_quantized = y_quantized(100:end);              % 去除初始不稳定值

%% 生成Verilog测试文件
fid = fopen('cic_input.txt', 'w');
fprintf(fid, '%x\n', y_quantized); % 十六进制格式写入
fclose(fid);
y_signed = double(y_quantized) - 8;  % 转换为有符号数（-8 ~ +7）

plot(y_signed);
grid on;

%% CIC积分器（96MHz时钟域）
integrator = zeros(N, length(y_quantized));
integrator(1, :) = mod(cumsum(y_quantized), 2^output_bit);       % 第一级积分
for i = 2:N
    integrator(i, :) = mod(cumsum(integrator(i-1, :)), 2^output_bit); % 多级积分
end

%% 降采样（每24点取1点）
downsampled = integrator(N, 1:D:end);

%% CIC梳状器（4MHz时钟域）
comb = zeros(N, length(downsampled));
comb_delay = zeros(N, length(downsampled));
comb(1, :) = mod((downsampled - [0, downsampled(1:end-1)]), 2^output_bit); % 第一级梳状
for i = 2:N
    comb(i, :) = mod((comb(i-1, :) - [0, comb(i-1, 1:end-1)]), 2^output_bit); % 多级梳状
end

%% CIC输出（最终结果）
golden_output = comb(N, :) - 2^(output_bit-1);

% 将golden_output转换为整数并写入文件
%golden_output_int = int32(golden_output); % 转换为32位整数
fid = fopen('golden_output.txt', 'w');
fprintf(fid, '%d\n', golden_output);
fclose(fid);

figure;
subplot(3,1,1);
plot(t, s, 'b', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Amplitude');
title('模拟波形');
xlim([0, 5e-6]);  % 显示前5μs细节
grid on;

%% 绘制CIC滤波后数字波形
subplot(3,1,2);
t_cic = (0:length(golden_output)-1) / (Fs/D); % 4MHz时间轴
%stem(t_cic, golden_output, 'r', 'MarkerSize', 4, 'LineWidth', 1.2);
plot(t_cic, golden_output);
xlabel('Time (s)');
ylabel('Amplitude (Digital)');
title('CIC滤波后输出（数字波形，4MHz）');
%xlim([0, 5e-6]);
grid on;

subplot(3,1,3);
plot(20*log10(abs(fft(golden_output))));

